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Abstract 

A  soft  ion  (SI)  model  of  1:1  electrolytes  (0.1  to  5  molar  concentration)  has  been 
studied  by  Grand  Canonical  Monte  Carlo  (GCMC)  simulations  and  hypemetted  chain 
(HNC)  integral  equation  theory.  Pair  correlation  functions,  osmotic  pressure,  and  mean 
activity  coefficients  have  been  calculated  as  a  function  of  the  repulsive  soft  core  potential 
r”*'  for  1/  =  9  and  12  (soft  ions)  and  u  =  <x>  (hard  ions).  Results  obtained  by  HNC  theory 
for  different  electrolyte  concentrations  and  soft  core  potentials  are  in  excellent  {^reement 
with  GCMC  simulations.  At  high  electrolyte  concentration  (  >  1.0  M),  the  thermodynamic 
properties  and  pair  correlation  functions  of  the  electrolytes  are  strongly  dependent  on  u. 
Activity  coefficients  obtained  with  a  soft  core  repulsion  potential  axe  in  good  agreement 
with  experimental  data  for  1:1  electrolytes.  Comparison  of  GCMC  results  obtained  using 
the  minimum  image  (MI)  and  Ewald  sununation  methods  indicate  that  the  less  expensive 
MI  method  gives  satisfactory  results  for  mediiun  to  high  electrolyte  concentration. 


I.  Introduction 


The  Debye-Hiickel  (DH)  [1]  theory  of  electrolytes  has  been  very  successful  in  describing 
the  structure  of  dilute  electrolytes.  The  finite  size  of  ions,  however,  limits  the  applicability 
of  the  DH  theory  to  extremely  low  ion  concentration.  One  of  the  most  often  tised  models  of 
electrolytes  which  takes  into  account  the  fibaite  size  of  ions  is  the  primitive  model  (PM).  In 
this  model,  the  ions  interact  through  a  repulsive  hard  sphere  core  potential  and  a  Coulomb 
potential.  Systematic  investigations  of  the  PM  electrolytes  using  approximate  theories  and 
simulations  indicate  that  an  oscillatory  structure  exists  in  the  ion  pair  correlation  functions 
at  high  concentrations  [3,5,9-13]  is  not  predicted  from  DH  theory. 

In  a  recent  report  [14],  Moleciilar  dynamics  (MD)  simulations  and  hypemetted  chain 
(HNC)  integral  equation  theory  were  used  to  study  the  effect  of  the  inverse  power  law 
intermolectdar  repulsion  l/r"*'  (with  v  —  4, 6, 9, 12, 24, and  oo)  on  the  structure  of  1:1 
electrolytes.  Two  model  systems  were  investigated.  The  first  (a)  is  a  modified  version  of 
the  (PM)  electrolyte  in  which  all  ions  were  soft  spheres  and  the  solvent  is  represented  by 
a  dielectric  constant  c.  This  model  is  referred  to  as  the  soft  ion  model  (SI).  In  the  second 
model  (b),  both  ions  and  solvent  molecules  were  represented  by  soft  spheres  eind  the  polar 
interactions  were  represented  by  a  dielectric  constant  c.  This  model  is  referred  to  as  the  soft 
ion  in  solution  (SIS)  model.  The  results  obtained  from  these  investigations  indicate  that 
the  repulsive  interaction  potential  plays  an  important  role  in  determining  the  structural 
and  dynamic  properties  of  electrolytes.  For  instance,  inclusion  of  the  solvent  imparted 
liquid-like  structure  to  the  ion-ion  pair  correlation  functions  and,  for  u  ^  oo,  gave  rise 
to  substantial  back  scattering  oscillations  in  the  velocity  autocorrelation  functions  with  a 
consequent  reduction  in  the  particle  self  diffusion  coeflScient.  Larger  values  of  v  reduced 
the  back  scattering  in  the  velocity  autocorrelation  functions.  The  dynamical  behavior  of 
the  ions  in  the  SI  model  is  similar  to  that  of  the  (PM)  electrolyte  {y  =  oo)  when  i/  >  9. 

In  the  studies  reported  here,  the  thermodynamic  properties  for  SI  model  of  1:1  elec¬ 
trolytes  are  obtmned  from  both  grand  canonical  Monte  Carlo  simiilations  (GCMC)  and 
HNC  theory.  Comparison  of  several  physical  properties  obtained  from  simulation  and 
theory  are  used  to  test  the  applicability  of  the  HNC  theory  over  a  range  of  electrolyte 
concentrations.  Ionic  activity  and  osmotic  coefficients  are  calculated  from  simulations  and 
therefore  chemical  potentials  can  be  obtained.  Comparisons  of  the  mean  ionic  iictivity 
coefficients  with  experimental  data  for  aqueous  solutions  demonstrate  the  validity  of  these 
soft  ion  model  potentials. 

II.  Model  systems  and  Interactions 

The  studies  are  carried  out  for  bulk  1:1  electrolytes  of  ion  concentrations  ranging  from 
0.1  to  5  molar.  The  unit  cell  for  the  GCMC  simulations  is  a  cubic  box  of  size  X  x  £  x  £  with 
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periodic  botmdary  conditions  in  all  three  dimensions.  The  interaction  potential  between 
ions  *  and  j  is  given  by  (14] 


<rij) 


6vd 


{d/rr  + 


iiSi 


(2.1) 


where  d  is  the  ion  size  pareuneter,  qt  is  the  ionic  charge  of  species  *  (value  tim^  the 
electronic  charge),  e  is  the  dielectric  constant  (taken  to  be  78.5),  and  Am  is  the  Madlung 
constant,  equals  1.74,  v  determines  the  range  of  repulsive  core  interaction.  This  repulsion 
comes  from  the  Pauli  exclusion  principle  of  the  filled  electronic  shells.  For  v  =  oo,  the 
potential  can  be  written  as 


rij  <  d 

Tij  >  d. 


(2.2) 


which  corresponds  to  primitive  model  of  electrolytes,  with  a  hard  core  of  diameter  equal 
to  d. 


III.  HNC  Theory 


In  a  bulk  electrolyte  solution,  the  structmal  properties  are  governed  by  the  Omstein- 
Zemike  (OZ)  equation 


hijiru) 


co(ri2)  + 

k 


J  Cikiri3)hkj{r32)dr3 


(3.1) 


where  n*  is  the  number  density  of  species  k{k  =  — ),  c<j(ri2)  and  ^i;(ri2)  =  gijiru)'- 1 

are  the  direct  and  total  correlation  functions.  In  the  hypemetted  chain  (HNC)  approxi¬ 
mation 

h{r)  =  c(r)  -j-  ln</(r)  -f  ^u{r)  (3.2) 

where  —  l/Ar^T,  ks  is  the  Boltzmann  constant  and  T  is  the  temperature.  The  OZ 
equation  becomes  [15] 

inyo(^i2)  =  / [/(ri3)y(ri3)  +  y(ri3)  -  1  -  Iny(ri3)]j*[(/»(r32)  +  l)y(r32)  -  l]*;dr3 

k  ^ 

(3.3) 

with 

y(r)  =  y(r)(/(r)  +  1)  (3.4) 

/(r)  =  —  1  is  the  Mayer  function.  Using  Newton  Raphson  method  with  an  initial 

guess  for  y(r),  eq.(3.3)  is  solved  with  300  nodes  and  node  spacing  of  0.02,  and  the  residue 
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The  excess  eaergy  and  osmotic  pressure  are  defined  as 


£^  =  2xnVViiX_,  /  Uij(r)gij(r)r'^dr  (3.5) 

>=i  Jo 

and 

PV/NkT  =  1  -  X"  E  E  r  (3.6) 

^  i=i  •'0 

where  u?(r)  =  r^,  Xi  is  the  mole  fraction  of  species  t,  and  n  is  the  number  density  with 
n  =  N/V. 

IV.  GCMC  Simulations 

The  Grand  canonical  ensemble  Monte  Carlo  (GCMC)  method  was  used  for  all  sim¬ 
ulations  reported  herein.  The  minimum  ims^e  method  (MI)  of  summing  the  Coulomb 
potentials  has  found  to  be  satisfactory  for  simulations  of  a  number  of  ions  of  ~  100.  The 
minimum  image  method  is  used  through  out  this  work,  unless  specifically  stated,  since  it 
is  considerably  less  expensive  than  the  Ewald  summation  method.  In  each  simulation,  the 
first  3  X 10^  steps  are  used  to  equilibrate  the  system  and  the  physical  properties  are  calcu¬ 
lated  from  the  following  equilibrium  simulations.  The  GCMC  simulations  are  performed 
as  following  [11,16]:  If  the  probability  of  accepting  a  trial  step  from  state  i  to  state  j  is 
fij,  then 

where  emd  Nf  sure  the  number  of  cations  and  anions  in  states  i  and  j  and 

B  =  2  [ln7±  +  In(nV’)  ]  ,  (4.2) 

with  7±  the  mean  ionic  sictivity.  The  Markov  chain  is  set  up  as  the  following 

fij  =  min(l,  fij/fji)  addition, 

fji  =  mtn(l,  fji/fij)  deletion,  (4.3) 

fij  =  mm[l,  exp(-y9l7ij)]  move. 

During  the  simulation,  the  probability  of  sulding  or  deleting  ions  is  varied  from  0.1- 
The  starting  msodmum  displacement  is  JD/4  smd  is  stdjiisted  in  the  first  10000  steps 
of  the  simulation  such  that  the  probability  of  accepting  a  psuticle  move  approaches  ~  50 
percent  if  possible.  Otherwise  the  maximum  displacement  is  the  value  which  is  fixed  at 
the  final  skdjustment.  The  number  of  particles  is  allowed  to  fluctuate  but  the  total  charge 
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in  the  simulation  is  kept  constant.  Therefore,  at  each  addition  or  deletion,  a  neutral  ion 
pair  is  added  or  deleted  to  maintain  charge  neutrality.  System  sizes,  the  number  of  steps 
executed  for  each  simulation,  and  simtdation  results  are  summarized  in  Tables  I-IV. 

For  the  primitive  model  simulations,  the  pressure  is  calculated  via  the  formula 


PV/NkT  =  ^  =  1  +  U/ZNkT  +  XiXjgij(d).  (4.4) 

»  j 


For  simulations  of  the  soft  ion  model,  which  has  a  continuous  repulsive  potential,  the 
pressure  is  calculated  using  the  virial  theorem 


^  =  1  + 


1 

3Nkf 


(4.5) 


V.  Results  and  Discussions 

Simvdation  and  HNC  integral  equation  theory  studies  of  three  different  core  potentials 
and  of  different  concentrations  will  be  presented  in  the  following.  To  quantitatively  test  the 
applicability  of  the  minimum  image  method  and  the  system  size  dependence  of  simulated 
thermodynamic  properties,  simulations  with  the  number  of  particles  from  20  to  216  with 
MI  method,  and  simulations  using  Ewald  summation  and  system  size  of  =  64  were 
performed  for  i/  =  9.  The  physical  properties  are  given  in  table  I,  and  the  pair  correlation 
functions  in  figures  1  and  2  for  C  =  1  molar  and  2  molar.  For  both  C  =  1  molar  and  C 
=  2  molar,  the  MI  method  gives  pair  correlation  functions  for  N  =  64  that  are  in  good 
agreement  with  results  obtained  from  Ewald  s\un.  The  dependence  of  the  pair  correlation 
frmctions  obtained  using  the  MI  method  on  the  system  size  is  also  shown  in  figtires  1  and 
2.  The  dependence  of  potential  energy  and  osmotic  pressure  upon  system  size  N  is  much 
weaker  than  that  of  the  ionic  activity  coefficient  as  demonstrated  in  Table  I.  If  N  =  64  is 
used  in  the  simulations,  the  MI  method  uses  considerably  less  computational  time  than 
the  Ewald  summation.  This  system  size  gives  a  good  compromise  between  accuracy  and 
cost. 

Figures  3(a-f)  is  the  calculated  and  simrilated  pair  correlation  functions  gij(r)  of  i/  =  9, 
for  electrolyte  concentrations  between  0.1  molar  to  5  molar.  These  figures  show  that 
the  HNC  theory  is  in  satisfactory  agreement  with  the  GCMC  simulations  within  this 
concentration  range.  As  the  concentration  increases,  the  difference  between  the  like  ion 
pair  correlation  functions  increases  slightly.  The  simulation  results  show  a  larger  charge 
inversion  at  higher  concentrations  in  the  second  layer  than  that  of  HNC.  At  C  =  4M,  there 
are  oscillations  evident  in  the  pair  correlation  functions  obtained  by  both  methods. 
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In  figvires  4(a-f),  one  can  see  even  better  agreement  between  the  simulated  and  calcu¬ 
lated  pair  distribution  functions  for  v  =  12.  Again  the  oscillations  in  the  pair  distribution 
functions  appear  at  high  concentrations  for  both  methods. 

Peiir  correlation  functions  of  the  PM  obtained  by  simulation  and  HNC  equation  are 
shown  in  figtires  5(a-f).  From  these  results,  it  is  clear  that,  the  HNC  integral  equation 
theory  works  better  for  larger  steeper  core  repulsion  (large  u).  To  compare  results  of 
the  PM  and  SI  models,  simulations  and  HNC  calculations  of  the  PM  electrolyte  using  an 
effective  heird  sphere  diameter  d  =  4.2  A  were  performed.  This  diameter  is  equal  to  the 
peak  position  of  in  the  SI  simulations  and  is  considered  to  the  effective  size  of  a  soft 
ion.  Comparisons  between  results  of  ions  which  have  the  same  effective  size  axe  believed 
to  be  more  meaningful.  Still,  the  rough  similarity  between  the  pair  distribution  functions 
obtained  from  the  two  models  exist  only  at  fairly  low  concentrations,  figures  6(a-d).  At 
higher  concentration,  the  PM  and  SI  models  give  very  different  results. 

Theoretical  investigations  show  that  the  Debye-Htickel  theory  breaks  down  at  /cd  ~  1 
[3-5,9-10]  for  the  PM,  where  k  is  the  inverse  Debye-Hiickel  (DH)  screening  length 

V  ribT 

At  ion  concentration  of  C  =  0.1  M,  the  DH  and  HNC  theories  give  almost  identical  results. 
At  C  =  5  M  the  DH  gives  an  unlike  ion  pair  correlation  function  peak  is  about  22%  less 
than  that  obtained  from  HNC  theory.  The  simulation  results  in  figures  5(e-f)  show  a 
charge  inversion  which  is  not  predicted  from  the  DH  theory.  Prom  Table  III,  we  conclude 
that  the  onset  of  the  oscillations  is  consistent  with  the  predictions  of  /cd  ~  1  [4]. 

The  thermodyneunic  properties  gives  a  much  clearer  picture  of  the  dependence  of  the 
SI  model  on  the  steepness  of  the  core  repulsion  piotential,  figures  7-9  and  Tables  I-IV.  At 
infinite  electrolyte  dilution,  all  three  model  potentials  {u  =  9, 12,  oo)  give  the  same  osmotic 
pressure  and  potential  energy  as  well  as  the  same  ionic  activity.  This  is  the  correct  ideal 
gas  limit.  At  finite  concentration,  the  thermodynamic  properties  are  a  strong  function  of  u. 
In  figure  9,  a  comparison  is  made  of  simulated  ionic  activity  coeflBccients  with  experimental 
values  of  several  1:1  electrolytes  [17-18].  The  mean  approach  of  ions  is  4.0  A  for  most 
of  these  salts.  Given  the  crude  description  of  the  core  repulsions,  the  soft  ion  model  gives 
very  good  description  of  the  ionic  activity  of  these  salts.  Using  the  PM  and  the  same  ionic 
size  parameter,  the  results  are  not  too  bad  for  KCl.  However,  the  results  obtained  from 
the  PM  seem  too  steep  to  characterize  most  of  the  salts. 
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Figure  Captions 

Figure  1.  Variation  of  the  pair  distribution  function  of  the  soft  ion  model  1:1  electrolyte  (t'  =  9) 
as  a  function  of  the  simulation  cell  sizes  (N  is  the  number  of  particles).  C  =  1  molar. 
Unless  specified,  the  simulation  is  done  by  using  minimum  image  (MI).  Keys  :  N  =  20, 
(o)  for  (+-),  (•)  for  (++);  N  =  64,  (A)  for  (+-),  (v)  for  (++);  N  =  216,  (+)  for 
(H — ),  X  for  (++);  N  =  64  using  Ewald  s\immation,  — )  for  (-^ — ),  ( - )  for  (++). 

Figure  2.  Variation  of  the  pair  distribution  function  of  the  soft  ion  model  1:1  electrolyte  (u  =  0) 
as  a  function  of  the  simulation  cell  sizes  {N  is  the  number  of  particles).  C  =  2  molar. 
Unless  specified,  the  simulation  is  done  by  using  minimum  image  (MI).  Keys  :  N  =  20, 
(o)  for  (+-),  (•)  for  (++);  N  =  64,  (A)  for  (+-),  (v)  for  (++);  N  =  216,  (+)  for 
(-1 — ),  X  for  (++);  N  =  64  using  Ewald  summation,  ( — )  for  (H — ),  ( - )  for  (++). 

Figures.  Comparison  of  pair  distribution  functions  for  the  soft  ion  model  1:1  electrolyte  for 
u  =  9  and  d  =  2.8  A  obtained  from  GCMC  simulations  and  those  from  HNC  integral 
equation  theory.  Keys  ;  GCMC,  (o)  for  (+-),  (•)  for  (++);  HNC,  (•  •  •)  for  (+-), 
( - )  for  (++). 

Figure  4.  Comparison  of  pair  distribution  functions  for  the  soft  ion  model  1:1  electrolyte  for 
t'  =  12  and  d  =  2.8  A  obtained  from  GCMC  simiilations  and  those  from  HNC  integral 
equation  theory.  Keys  :  GCMC,  (o)  for  (H — ),  (•)  for  (++);  HNC,  (•  •  •)  for  (H — ), 
( - )  for  (++). 

Figure  5.  Comparison  of  pair  distribution  functions  for  the  PM  model  1:1  electrolyte  from 
GCMC  simulations  and  the  HNC  integral  equation  theory  (d  =  2.8  A  ).  Keys  ; 
GCMC,  (o)  for  (+-),  (•)  for  (++);  HNC,  (•  •  -)  for  (+-),  ( - )  for  (++). 

Figure  6.  Comparison  of  pair  distribution  function  of  the  soft  ion  model  1:1  electrolyte  {  u  =  9, 
d  =  2.8,  and  u  =  oo,  d  =  4.2)  obtained  from  GCMC  simulations  with  HNC  results 
for  the  primitive  model  1:1  electrolyte  with  an  effective  diameter  d  =  4.2  X  (v  =  oo). 
Keys  :  GCMC,v  =  oo,d  =  4.2,  (o)  for  (4 — ),  (•)  for  (++);  GCMC,u  =  9,d  =  2.8, 

( - )  for  (+_),  (_)  for  (++).  HNC,  u  =  oo,d  =  4.2,  ( - )  for  (+-),  (•  •  •) 

for  (++). 
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Figure  7.  Osmotic  coefficients  ^  for  a  1:1  electrolyte  with  i/  =  9, 12,  and  oo  as  a  function  of  the 
concentration.  Lines  for  GCMC,  points  for  HNC. 

Figvire  8.  Potential  energy  U  for  a  1:1  electrolyte  with  v  =  9, 12,  and  oo  as  a  function  of  the 
concentration.  Lines  for  GCMC,  points  for  HNC. 

Figure  9.  Comparison  of  ionic  activity  coefficient  (ln7±)  for  1:1  electrolyte  (i/  =  9, 12,  and  oo) 
as  a  function  of  the  concentration  with  experimental  results.  The  curve  corresponding 
to  df.  is  from  simulations  using  an  effective  ionic  radius  d*  =  4.2  A  for  i/  =  9. 
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Table  I.  Parameters  and  results  of  the  HNC  calculations  and  GCMC  simulations  of  a 
l*.l  electrolyte  at  different  concentrations,  u  =  9.  The  side  length  of  the  cubic  simulation 
box  is  i,  the  average  number  of  ions  is  <  N  >,  U  is  the  average  potential  energy,  C  is 
the  bulk  ionic  concentration  in  molar,  iV,  is  the  number  of  steps  run  in  the  simulation,  <i> 
is  the  osmotic  coefficient,  and  7±  is  the  ionic  activity,  B  is  the  input  parameter. 


~wr 

B 

<N> 

-U/NkT 
MC  HNC 

4 

MC 

HNC 

C{M) 

ln7± 

Ns 

10® 

81.0 

6.410 

63.87 

0.278 

0.274 

0.933 

0.935 

0.100 

-0.259 

5 

81.00® 

6.413 

64.23 

0.278 

0.933 

0.100 

-0.263 

20 

47.37 

6.203 

63.24 

0.453 

0.455 

0.947 

0.948 

0.494 

-0.352 

5 

47.37® 

6.227 

64.10 

0.454 

0.948 

0.501 

-0.353 

20 

25.51 

3.978 

20.59 

0.549 

0.536 

1.004 

1.007 

1.029 

-0.3424 

5 

29.84 

4.905 

32.18 

0.540 

1.003 

1.006 

-0.3257 

5 

37.60 

6.345 

65.30 

0.534 

1.008 

1.020 

-0.3131 

2 

47.37 

7.678 

127.38 

0.527 

1.006 

0.995 

-0.3152 

1 

56.39 

8.744 

215.09 

0.525 

1.006 

0.996 

-0.3058 

1 

37.60® 

6.279 

63.44 

0.531 

1.003 

0.991 

-0.317 

20 

20.25 

4.345 

20.38 

0.619 

0.595 

1.154 

1.158 

2.038 

-0.1489 

10 

23.68 

5.301 

32.39 

0.603 

1.153 

2.024 

-0.1342 

3 

29.84 

6.691 

64.22 

0.592 

1.154 

2.007 

-0.1232 

2 

37.60 

8.078 

125.51 

0.585 

1.149 

1.961 

-0.1002 

1 

44.76 

9.124 

213.73 

0.582 

1.154 

1.979 

-0.1094 

1 

29.84® 

6.671 

63.49 

0.592 

1.151 

1.984 

-0.122 

20 

29.84 

9.350 

128.37 

0.578 

0.586 

1.554 

1.564 

4.012 

0.5133 

2 

23.68® 

7.958 

64.42 

0.586 

1.552 

4.026 

0.507 

20 

29.84 

10.624 

160.98 

0.536 

0.545 

1.807 

1.823 

5.031 

0.9239 

3 

*  Ewald  sum  method  is  used. 
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Table  II.  Same  as  in  Table  I  with  t/  =  12 


i(A) 

B 

<N> 

-U/NkT 
MC  HNC 

<i> 

MC 

HNC 

C(M) 

In7± 

Ns 

10® 

81.0 

6.410 

65.83 

0.292 

0.289 

0.926 

0.927 

0.103 

-0.2889 

2 

47.37 

6.118 

63.81 

0.486 

0.489 

0.922 

0.923 

0.498 

-0.4035 

2 

37.60 

6.142 

65.41 

0.584 

0.591 

0.957 

0.959 

1.022 

-0.4167 

2 

29.84 

6.288 

63.63 

0.674 

0.681 

1.049 

1.054 

1.988 

-0.3162 

2 

29.84 

8.379 

126.56 

0.734 

0.748 

1.304 

1.309 

3.955 

0.0419 

2 

29.84 

9.421 

162.94 

0.740 

0.754 

1.481 

1.493 

5.092 

0.3103 

2 

Table  III.  Same  as  in  Table  I  with  i/  =  oo  and  d  =  2.8  A. 


1(A) 

B 

<N> 

-U/NkT 
MC  HNC 

</> 

MC 

HNC 

C(M) 

-H 

1 

Ns 

10* 

nd 

81.0 

6.351 

64.60 

0.313 

0.307 

0.912 

0.918 

0.101 

0.299 

3 

0.2919 

47.37 

5.931 

63.68 

0.547 

0.542 

0.877 

0.889 

0.497 

0.495 

3 

0.6475 

37.60 

5.835 

64.93 

0.676 

0.671 

0.881 

0.896 

1.015 

0.563 

3 

0.9253 

29.84 

5.771 

64.90 

0.815 

0.809 

0.913 

0.936 

2.028 

0.594 

3 

1.3079 

29.84 

7.366 

128.27 

0.958 

0.956 

1.019 

1.050 

4.008 

0.478 

4 

1.8387 

29.84 

7.971 

15.11 

1.006 

1.004 

1.078 

1.115 

4.941 

0.384 

4 

2.0413 

Table  IV.  Same  as  in  Table  I  with  u  =  oo  and  d  =  4.2  A. 


X(A) 

B 

<N> 

-U/NkT 
MC  HNC 

MC 

HNC 

C{M) 

b7± 

N, 

10® 

81.0 

6.451 

63.92 

0.274 

0.270 

0.942 

0.944 

0.100 

-0.2388 

2 

47.37 

6.415 

65.29 

0.464 

0.458 

0.979 

0.989 

0.510 

-0.2783 

2 

37.60 

6.611 

64.33 

0.560 

0.553 

1.079 

1.084 

1.005 

-0.1651 

2 

32.14 

7.384 

64.59 

0.671 

0.659 

1.320 

1.345 

2.019 

0.2172 

2 

11 


0.0  2.0  4.0  6.0 

Figure  1  Distance  r/d 


Pair  Correlation  Functions  g(r) 


0.0  2.0  4,0  6.0 

Figure  2  Distance  r/d 


Pair  Correlation  Function  g(r) 


12  3  4  1  2  3  4 

Figure  3  Radial  Separation  r  (A) 
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